Viscoelasticity of two-layer-vesicles in solution 
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The dynamic shape relaxation of the two-layer-vesicle is calculated. In additional to the undu- 
lation relaxation where the two bilayers move in the same direction, the squeezing mode appears 
when the gap between the two bilayers is small. At large gap, the inner vesicle relaxes much faster, 
whereas the slow mode is mainly due to the outer layer relaxation. We have calculated the vis- 
coelasticity of the dilute two-layer-vesicle suspension. It is found that for small gap, the applied 
shear drives the undulation mode strongly while the slow squeezing mode is not much excited. In 
this limit the complex viscosity is dominated by the fast mode contribution. On the other hand, 
the slow mode is strongly driven by shear for larger gap. We have determined the crossover gap 
which depends on the interaction between the two bilayers. For a series of samples where the gap 
is changed systematically, it is possible to observe the two amplitude switchings. 



I. INTRODUCTION 

Vesicle dynamics has long been investigated both ex- 
perimentally and theoretically [I}j6]. The decay rate of 
the thermally excited shape fluctuation provides infor- 
mation of vesicle elasticity. The knowledge of the micro- 
scopic relaxation can also be used to predict the macro- 
scopic rheological property of the vesicle solution [6]. 

The vesicles are mostly non-equilibrium system. The 
external perturbation, be it thermal, electrical, sonica- 
tion [7J IB] , or flow [HI HH] ; transforms the lamellar struc- 
ture into the vesicles. In such situation, one often ob- 
tains mixture of uni-lamellar vesicle and multi-lamellar 
vesicles (MLV) of various sizes. However, in sharp con- 
trast with the very detailed calculations and scattering 
experiments on the unilamellar vesicles, there are rela- 
tively few studies of the MLV dynamics in the literature. 
In this paper, we consider the two-bilayer vesicle which is 
the simplest MLV. We calculate the dynamic relaxation 
rates, as well as the rheological response of the two-layer- 
vesicles. From the experimental point of view, it is hard 
to prepare vesicles with exactly two bilayers. Nonethe- 
less, such a concrete calculation provides a clear picture 
of how the interaction between the two membranes affect 
the relaxation rates, as well as the squeezing (lubrication) 
flow which arises exclusively in MLV. Since our method 
can be extend to vesicles with more than two bilayers, a 
more specific calculation can be performed similarly. 

The squeezing dynamics of the sandwiched solvent be- 
tween the two layers have been analyzed in two related 
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problems. This relaxation process of the flat lamellar 
system was calculated by Brochard and de Gennes as 
the "slip mode" [TT], or later measured and called as 
"baroclinic mode" [P21 - fLl] . In the soap film system, the 
"squeezing mode" dispersion has been calculated and 
measured |15[ 116] . In this work, we obtain the similar 
relaxation in which the solvent is squeezed between the 
two adjacent bilayers to relax the bilayer curvature en- 
ergy and the mutual interaction energy between the two 
bilayers. As this mode arises only in MLV, we are inter- 
ested in its dispersion relation and its coupling strength 
with the applied shear. Due to the strong lubrication 
resistance, the squeezing mode is often the slowest re- 
laxation mode. This makes the squeezing mode an im- 
portant candidate for the vesicle rheology. However, not 
every relaxation mode is equally excited by the applied 
shear. Therefore whether the shear can drive the squeez- 
ing mode with a large amplitude is an important problem 
to consider. In this work, we will try to build our under- 
standing of the coupling strength through the concrete 
calculation. 

Based on the analysis below, we find that the squeezing 
mode makes the dominant contribution to the complex 
viscosity for strongly interacting bilayers. Interestingly, 
when the gap between the two bilayers is smaller than 
a characteristic crossover gap, the fast undulation mode 
becomes the dominant mode for the complex viscosity. 
We find that the crossover gap gets very small when the 
bilayer interaction is strong. On the opposite limit where 
the bilayer has extremely weak interaction, the crossover 
gap becomes comparable to the inner vesicle radius. 

Below in Sec.|Hj we define our model. Section [Hi] dis- 
cusses the elastic force. In Sec. [W] and Appendices, the 
flow resistance between the bilayers is solved. In Sec. |VJ 
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we shall analyze the relaxation rate. In Sec. VI we con- 
sider the viscoelasticity of the dilute vesicle suspension. 
Finally in Sec. |VII| we summarize and discuss our results. 



II. THE MODEL 




We consider a two-layer-vesicle with two bilayers lo- 
cated at the spheres with mean radius r\ and r-i as shown 
in Fig. [I] The deformation of the bilayer shapes are de- 
scribed by the radial layer displacements u\ and u-i rela- 
tive to the two reference spheres respectively. Both of the 
bilayer membranes are surrounded by a solvent of viscos- 
ity r\. Let a denote the area per molecule projected on 
the reference sphere, and ao the averaged projected area 
per molecule. For each membrane, we define the dimen- 
sionlcss surface density <f> n (n = 1, 2) given by the ratio 
ao/ a. Notice that <p n becomes unity at equilibrium. In a 
fluctuating vesicle, </>„ is not uniform in general. In this 
work we propose the free energy of a two-layer-vesicle as 
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where B is the layer compression modulus, j n the surface 
tension, k the bending modulus, and E the area stretch- 
ing/compression modulus. Here both n and E are taken 
to be the same for the two membranes. Moreover, dfi is 
the differential solid angle, and the surface area element 
is approximately given by dA n m [1 + (Vj_u„) 2 /2]r 2 dO. 
The mean curvature H n is given by 
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(2) 



up to linear order in u n . 

Two comments should be made about the first term of 
Eq. 0. The compression modulus B should be a func- 
tion of r\ and In principle a microscopic statistical 
model for MLV should provide the functional form. Here 
we focus on the layer dynamics, therefore to build such 
a model is beyond the scope of this work. As a rough 



estimate for discussion, below in subsection III Bl we will 



use the B, which is calculated from the flat layers, and 
replace the distance d between two layers by r 2 — T\ . This 
approximation is justified for d <C 7~2, while it may devi- 
ate considerably when d is comparable to r 2 . 

The interaction term proposed here is proportional to 
(u2 — ui) 2 , which arises naturally at small d from the 
square of the strain. At large d, the curvatures and 
the areas for the two bilayers are very different. Pre- 
sumably a microscopic model may derive a more suit- 
able weighted interaction energy, which is proportional 
to [u 2 — {r\/r2)^Ui\ 2 with a weighting exponent f3. Below 
we present the calculation without this extra weighting 



FIG. 1. The vesicle consists of two bilayers with the radius 
ri and T2- The surrounding solvent viscosity is rj everywhere. 



factor. Nonetheless, the calculation procedure is exactly 
the same for the case (3 ^ 0. 

At the dilute phase boundary of the lamellar phase, 
the surface tension j n should vanish. Inside the lamellar 
phase, the tension depends on the applied osmotic pres- 
sure. To simplify the discussion, we shall only consider 
the case where the surface tension vanishes 7„ = 0. In 
fact, the calculation with finite surface tension can be 
carried out in the same way. The layer displacement u n 
is related to the radial velocity v r at r = r n by 



du n 
~dt 



v r (r n ), 



(3) 



where we neglect the solvent permeation. The surface 
density obeys 



d(f> n _ 2v r (r n ) ^ 
dt r n 



-V_L-[vi(r n )0 n ], (4) 



where Vj^ is the two-dimensional (2D) surface derivative 
and v±(r n ) is the tangential velocity at r = r n . Here we 
have neglected the amphiphile exchange flux between the 
neighboring bilayer (the amphiphile permeation). 

The flow field obeys the Stokes equation which is pre- 
sented here as the form 



7?V 2 v 



Vp = 0, 



(5) 



where 77 is the viscosity, and p the pressure. Here we 
write down the force balance conditions that are satisfied 
on the two layers. The normal force balance is given by 



6F 

6u n 



+ <*rr(r+) -<T rr (r n ) = 0, 



(6) 
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where a rr = —p + 1r\d r v r [r being the radial distance), 
and the superscripts + and — indicate that the stresses 
are evaluated at the exterior and interior of the bilayer, 
respectively. On the other hand, the tangential force bal- 
ance is given by 



SF 



+ <rx{r+) -<rx(r n ) = 0, 



(7) 



where x„ is the 2D tangent displacement of the layers, 
and the 2D stress is defined as a± = 9a r g + (pa rip (9 and 
ip being the polar and azimuthal angles, respectively). 

The viscoelastic response of a dilute two-layer-vesicle 
suspension can be calculated by considering the stress 
response to an external flow at large distances from the 
vesicles 



v°°(r,t)=rV [r 2 Y 20 (9^)] e 1 ' 



(8) 



where T is the strength of the elongational flow, Y[ m (9, <p) 
are the spherical harmonics, and uj is the angular fre- 
quency. Each suspended vesicle contributes to the aver- 
aged stress. In the dilute solution, the effective complex 
viscosity is calculated as [5J IT?] 



\r\Yr\ ' 



(9) 



where <f> v = (4n/3)cr 2 is the volume fraction occupied 
by the vesicles (c is the number density of the vesicles) , 
and p\\ ' s the frequency dependent complex coefficient 
of the spherical harmonic expansion of the pressure (see 
Eqs. (A10) and (|B6|) later). The effective viscosity can 



also be expressed in terms of the complex modulus as 
G*=iur)*. 



III. THE ELASTIC FORCES 

A. Force expressions 

The elastic forces are calculated by evaluating the 
derivatives —8F/8u n and —SF/Sx n . From Eq. Q, we 
see that the bilayer tangential displacement and the ra- 
dial displacement produce the first order surface density 
perturbation as 



2u n 



The tangential elastic force is given by 
5F 



(10) 



(11) 



where the tension perturbation Sj n is induced by the 
surface density perturbation, i.e., Sj n = —E{(j) n — 1). 
Notice that the 2D gradient has the usual component 
form; 



V_l(<57„) = 9 — — + V^^-p.—^. • 

r c)9 r sin t) Oip 



(12) 



Up to linear order in u n and Sj n , the normal force on 
bilayer 1 is given by 
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whereas that for bilayer 2 is 
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In the large stretching modulus limit (E — > 00), the ten- 
sion perturbations Sjn become Lagrange multipliers, so 
that they ensure that the right hand side of Eq. ( |To| van- 
ishes. By using Eqs. Q and (11), the values of Sj n are 
determined from the viscous stress a± (see Appendix C). 



B. The bilayer interactions 

There are several interactions which contribute to the 
layer compression modulus B. In this subsection, we 
indicate their physical origins and give the simplest for- 
mulae to describe them. For charged bilayers, the elec- 
trostatic interaction, together with the counter-ion and 
co-ion entropy produce the free energy per area 



64C s fc B T + 2 f q$ 
tanh ' 



(15) 



where C s is the salt concentration, qip the potential en- 
ergy of the counter-ion q at the surface, kd the inverse 
of the Debye screening length, d = r 2 — r\ the distance 
between the two surfaces, and k^T the thermal energy. 

The van der Waals attraction potential per unit area 
can be calculated by summing the dipoles to obtain 
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where A is the Hamaker constant, and 5 is the membrane 
thickness. A more complicated Lifschitz theory calcu- 
lation can provide a more accurate description [TSHU]- 
The sum of these two interactions consists the standard 
DLVO theory. When the electrostatic repulsion and van 
der Waals attraction stabilize MLV, B can be evaluated 
from 



B 



d 2 (V c + V vdW ) 
' dd 2 



(17) 



For a flexible bilayer where the undulation entropy 
depends strongly on the membrane separation d, Hel- 
frich estimated the free energy per area and B by a self- 
consistent argument to obtain [THl H2] 
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where cq is an order unity numerical constant. In the 
following examples below, we use Co = 36 /tt 2 . When the 
van der Waals interaction becomes important, a more 
elaborate calculation is required to combine the undula- 
tion entropy and the van der Waals attraction [53]. 



IV. THE NORMAL FORCE BALANCE 

Using the spherical harmonic expansion, the Stokes 
equation can be solved in terms of the radial functions 
which arc simple polynomial of the radius. The detailed 
calculation is presented in the Appendix A. In general, 
the solution depends on the boundary velocities, which 
are the velocities of the membranes and the applied ex- 
ternal shear. In the usual case where the bilayers have 
a large stretching modulus E, we can simplify the dis- 
cussion by considering the large stretching modulus limit 
E — > oo. In this limit, the surface density approaches 
unity, and the combinations E(<p n — 1) become the La- 
grange multipliers to ensure that surface densities are 
constant. Hence the constant E can be eliminated. The 
detailed calculation is presented in Appendix B. At the 
boundaries, this limit also turns the tangential velocity 
into the function of the normal velocity. This means that 
the normal stresses cr rr (r+) and a rr (r~), as well as the 
tension perturbation 8^f n are all proportional to the bi- 
layer normal velocities and the external shear. The de- 
tailed calculation is presented in Appendix C. 

The bilayer displacements and velocities are expanded 
as the sum of the spherical harmonics Yi m (9,tp) times 
the time varying amplitudes, i.e., u n (9, tp) = u n Yi m (9, ip) 
and v r (r n , 9, tp) — v n Yi m (9 , tp) . Without cluttering the 
notation with the angular quantum numbers I and m, 
hereafter we use u n and v n to express the time varying 
amplitudes for an arbitrary set of (l,m). After perform- 
ing some calculations to express a rr and Sj n by v n (see 



Eqs. |C1J), |C5J), |B4], (|B6|) and (|B9j)), we find that the 
normal force balance Eq. ([6| at r% and r 2 can be written 



r 2 u 2 



D 



20r 1 Te^ t 5 l2 5 mO e 2 , (19) 



where e 2 = (0,1). Here we prefer to use the variables 
r\u\ and r 2 u 2 which make the matrices E and D sym- 
metric. With this variable choice, the free energy per 
solid angle is then the quadratic form of E. The compo- 
nents of the matrix E are given by 



E n = 



B(l-p 3 ) ( 7l r 2 2 p 2 + kL 2 )(L 2 - 2) 



3r 2 3 (l - p) V 
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where have defined p = i\/r 2 (< 1) and 

f 2 = 1 9 ( ■ f) d \ 1 92 
sin9d9 \ 39 ) sin 2 9 dp 2 ' 

Whereas the components of the matrix D are 

T^ + ^+l) 



(21) 



G (! 2 + l) 
[-(l + l) 2 (4l 2 -l)p 2l+1 

+ (l 2 -l)(2l-l)(2l + 3)p 21 - 1 
+ (l-l) 2 (4l 2 + 8l + 3)p 2l ~ 3 
+ (8l 2 + 81- 4) P - 2 ] , 



G (l 2 + l) 
[(2l 3 + 5l 2 + l-2)(p 3l+2 -p 1 - 1 ) 
+ (2l 3 + l 2 -3l)(p l+1 -p 31 )] , 
D 2 i = D 12 , 

_ V 4 l+1 (2l + l) 
° 22 - G Q (l 2 + l) X 

[-(l + 2) 2 (4l 2 -l)p 2l + 4 
+ 2l(l + 2)(2l~ \){2l + 3)p 2l+2 
+ l 2 (Al 2 + 8l + 3)p 21 

+ (8l 2 + 81- 4)p] , (22) 

G = [4p + 4p 4/ + 3 - (21 + l) 2 p 2l+i 

-(6 - 8^ - 8l 2 )p 2l+2 - (21 + l) 2 p 21 ] . (23) 

Notice that the components of D are the product of r//r 2 
and the functions of p and I. 

When the two bilayers are well separated, i.e., n <C r 2 , 
they are not hydrodynamically coupled. In this limit, 
D\ 2 and D 2 \ become small, and we recover the isolated 
vesicle damping given by [1HB] 



where 



r)(2l + l)(2l 2 + 21-1) 
r 3 J(l + l) 



(24) 



For the more general bilayer interaction which is pro- 
portional to B[u 2 — (ri/r 2 )^u{\ 2 , the similar calculation 
will give rise to an extra factor p 2 @ for the B term in En, 
and an extra factor p@ for the B terms in E\ 2 and -E 2 i. 
The form of -E 22 is not affected. 



V. RELAXATION SPECTRUM 

The relaxation rates, denoted as ftj (j — 1,2), are 
the eigenvalues of the matrix D _1 • E. The normal force 
balance condition Eq. Q on the two bilayers gives the 
eigenvector equation 
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where dj- are the (non-orthogonal) right-eigenvectors. To 
work with the orthogonal eigenvectors, here we define the 
symmetric decay rate matrix 



M = T>- 1/2 -E-D- 1/2 , 
which has the (normalized) eigenvectors dj 



M d, 



(26) 



(27) 



Here the eigenvalues are the same as in Eq. ( 25 ) , and d^ 



is proportional to D 1 / 2 d?'. We can decompose the decay 



rate matrix as 



M = ftididi + ft 2 d 2 d 2 , 



(28) 



which will be used in the later discussion. Hereafter the 
faster and the slower rates are denoted by f?i and fi 2 , 
respectively. 

The full expression of the decay rates are straightfor- 
ward, but complicated. Hence we shall obtain some sim- 
plified expressions to gain some understanding on the na- 
ture of the relaxation. We first discuss the fast mode fii. 
When p = 1, the fast undulation mode has a well-known 
dispersion relation [21 IH [6] 



fin 



2n(l-l)l 2 (l + l) 2 (l + 2) 
777-3(2/ + 1)(2Z 2 + 2Z — 1) ' 



(29) 



where the bending modulus is doubled as there are two 
bilayers. For p slightly less than unity, a series expansion 
can be made if desired. For p <C 1, the decay rate can be 
approximated by 



En 



(30) 



This corresponds to a simple picture that the fast mode 
consists mainly the inner bilayer relaxation, whereas the 
outer bilayer does not move much. 

Next we consider the slow mode fi 2 . When p m 1, the 
slow mode can be approximately obtained by the series 
expansion of 1 — p. Based on the relation fii + fi 2 = 
tr(D -1 • E), in which 2 only starts from the second 
order term, we get f2i ~ tr(D _1 • E) for the zeroth and 
the first order terms. Because the product fiif2 2 is given 
by det(D _1 -E), we use the leading two terms of det(D • 
E) /f2i to obtain 



Bltl + l), 

\ ' (l-p) 2 (2-p) 

+ K(l-l)l 2 (l + l) 2 (l + 2) _ 3 
24777-2 



(31) 



The leading term (B/12rf)(l - p) 2 l(l + 1) is similar to 
the slip mode of the planar smectic, which has the decay 
rate {B /12rf)d 2 q\ : where q± is the wave vector projected 
on the bilayer plane This expression is also analo- 
gous to the squeezing mode dispersion obtained for soap 



film [T51 [TB] . In the other extreme of p 
mode has the dispersion such that 



0, the slow 



[B + 3(1- 1)1(1 + l)(l + 2))K/rl] 1(1 + 1) 
2 ~ 3r?(2f + 1)(2Z 2 + 2Z- 1) ' [ ' 

If we set the shear perturbation as I = 2, the numerator 
contains a factor B + 72k/ Y 2 . Hence the dimensionless 
parameter Br\jK becomes important when its value is 
much larger than 72. 

We now discuss another approximate expression for 
the slow mode valid for the whole range of p. Since the 
sum of the two decay rates fli + fi 2 = tr(D _1 ■ E) is 
dominated by the fast mode for any p, the slow mode 
can be approximated by the ratio 



n 2 



dct(D- 1 • E) 



tr(D" J • E) 



E11E22 — EyiEi\ 



E\\D 2 2 — E 12 D 2 i + E 22 Dii — E 2 \Di2 ^ 

When the inner bilayer 1 relaxes fast, one can apply the 
adiabatic approximation to the fast relaxing inner mem- 
brane. The approximated expression has even a simpler 
denominator given by 



E22 — E2i(Ei2/En) 



(34) 



D22 — D2i(Ei2/ En) 

In Fig. [2j we plot the decay rates as a function of p = 
ri/r 2 for 1 — 2, keeping a constant B. The two solid 
lines represent the numerical calculated decay rates. The 
upper one is the fast mode fii which has the limit fi at 
p = 1. It coincides with the approximation Eq. (30) 



(dotted line) for p < 0.6. The lower solid line represents 
the slow mode fi 2 . The approximate slow rate Eq. (34) 



follows qualitatively the exact value of fi 2 for the full 
range of p. It also coincides with Eq. (32) in the limit 
of p = 0. Incidentally the inner layer relaxation Eq. (30) 
also fits the slow mode CI2 at p w 1. 

Fig. [2] is useful to illustrate the nature of the relax- 
ation for the fast and slow modes. In a real system, B is 
a function of r\ and r 2 (or p and r 2 ), which is hard to be 
kept as a constant. One will make a plot with a function 
for B which is suitable for the specific MLV system. If 
one uses the flat layer formulas described in subsection 
III B as approximations, caution should be kept for the 
accuracy at small p. Given an accurate function for B, 
together with the proper weighting factor (to combine 
p 2 P and pP factors in E), the above decay rate approxi- 
mations should work better. 



VI. VISCOELASTICITY 

We now consider two-layer- vesicles under the external 
oscillatory shear flow. Rearranging Eq. ([9]) together with 
Eq. (B6), we obtain 



v* hi ~ 1 



V2 



r 2 re l 



(35) 
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FIG. 2. The scaled relaxation rate Qj/Qo as a function of 
the dimensionless size ratio p — ri/r% between the two layers. 
Here Qo is the vesicle relaxation rate with the rigidity 2k given 
by Eq. ( |29[ ). The two solid lines represent the two vesicle 
relaxation rates fii and Q.2 obtained numerically. The dotted 
line is an approximation for the fast mode given by Eq. ( 30 1. 



The dashed line repr esen ts the approximate slow relaxation 
rate ^2 given by Eq. (34 1. In this plot, we set B — 600 J/m 3 , 
k = k B T = 4 x 10~ 21 J, 77 = 10~ 3 Pa-s, r 2 = 10"" 
so that Brljn = 150. 



I = 2, 



FIG. 3. The shear coupling strength S defined by Eq. (371 as 
a function of the dimensionless size ratio p between the two 
layers. 



It is apparent that the dilute hard sphere limit is recov- 
ered when V2 — 0. By putting u n — v n /iuj and solving 



Eq. (19) for V2, we obtain 



5 20r/ A 

2 " iW 7F e2 



(E- 



■jD) 



e 2 , 



(36) 



where one should set / = 2 for the matrixes E and D. 
Here £2 appears twice, because the shear directly affects 
the bilayer 2 through Eq. (19 1, and the velocity of the 



bilayer 2 carries the stress contribution of the vesicle ac- 
cording to Eqs. ([9]) and (B6). 

Considering the high-frequency limit, we now define 
the dimensionless shear coupling strength S and the 
shear deformation unit vector s by 



\ r 2 



1/2 



D-i/a 



e 2 - 



(37) 



Notice that S and s only depend on the ratio p as long as 
we fix to I — 2. As shown in Fig. [3} S is weakly depended 
on p because its value varies only between 2.1 and 2.19. 
Hence the right hand side of Eq. ( [36] ) varies between 0.31 
and 0.4, whereas the hard sphere case gives 2.5. This 
means that, at high frequencies, the coupling of a vesicle 
to the flow is weaker compared to that of a hard sphere. 
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P P 



FIG. 4. The relative viscosity strengths defined by Eq. (391 



as a function of the dimensionless size ratio p between the 
two layers. The solid and the dashed lines corresponds to 7/1 
and r)2, respectively. We chose the value Br\l k = 150. The 
crossing of the two modes occurs at p* « 0.733. 



FIG. 5. The polar angle a (divided by 7r) of various 2D vectors 
as a function of the dimensionless size ratio p between the two 
layers. The solid line represents s, while the dotted and the 
dashed lines represent di and &2, respectively. We chose the 
value Br\ / K = 150. 



The high frequency viscosity asymptotically ap- 
proaches 



00 

7] = 7] 



1 



s 



rf -r)S<p v , 



(38) 



where rf = r][l + (5/2)0 v ]. Using the eigenmode decom- 
position I = didi + d2d2, we can separate the viscosity 
contributions from each mode as 



so that we have 



T]S = 1] 1 + f] 2 . 



(39) 



(40) 



Then the full expression for the complex viscosity is given 
by 



rf (ui) = if - (f) v ^2 Vj 
3=1 



iu> + CI a 



(41) 



The relative strength of the two modes are shown in 
Fig. [4j Alternatively (but equivalently) , the complex 



modulus can be expressed as 



2 

G*(uj) = iuir) 00 + (j) v y" G 3 U 



3=1 



■01} 



(42) 



where the two modes have the positive amplitudes Gj — 

Even for different bilayer interaction strength B, we 
always find that the slow mode has the larger viscosity 
amplitude at p s=s 0, and smaller one at p ~ 1. The latter 
is reasonable because the shear perturb the two bilayers 
similarly for p « 1, where s is along the (1, Indirection. 
In terms of the polar angle a between s and the first axes 
on the (rfui, r|it2)-plane, the (1, Indirection corresponds 
to a — 7r/4. Notice that this direction also corresponds 
to the undulation eigenmode direction. Since the sum 
of the two viscosity amplitudes is roughly a constant as 
shown in Eq. ( 40 ) , the slow squeezing mode takes a small 
viscosity amplitude 772 for p w 1. In the other limit of 
p « 0, the shear mainly perturbs the outer layer. In this 
limit, we have § = (0, 1) or a = n/2, and the slow mode is 
due to the outer layer relaxation. As a result, rj 2 becomes 
the dominant contribution for psiO. In Fig. [5j we plot 



VII. SUMMARY AND DISCUSSION 



0.9- 



0.8- 



^0.7 



0.6- 
0.5 
0.4 



fast mode 



slow mode 



10 



10 



10 

Br 2 3 /K 



10 



FIG. 6. The crossover size ratio p* as a function of the di- 
mensionless bilayer interaction Br\jK. For p > p* , the fast 
mode is the dominant contribution to the viscosity, while the 
slower relaxation is dominant for p < p* . 



the angle a as a function of p. The shear vector s always 
coincides with the fast mode at p — 1 and the slow mode 
at p = 0. Therefore the mode switching always takes 
place between < p < 1. It is also worth mentioning that 
Figs. Q and (J5J), like Fig. ([2]), are plotted by assuming a 
constant B. These plots are used to analyze the behavior 
of the mode amplitudes. The more physical plots will be 
the ones using a suitable function for B. 

Here we define a crossover size ratio p* at which the two 
modes have the same viscosity amplitude, i.e., r\\ = ij 2 . 
In Fig.[6j we plot the calculated p* as a function of Br\j n. 
When Br 2 /n <C 1, the crossover happens at p* ss 0.52. 
This means that for non-interacting case B = 0, the 
crossover happens when d = r 2 —r\ is roughly the same as 
n. When Bt 2 /k ^> 1, on the other hand, p* approaches 



unity. As mentioned in Eq. (32 1, the parameter Br 2 /n 



will have a noticeable effect only when it is greater than 
72. At the lower right corner of Fig. [6j where the layer 
interaction is strong and the layer separation is not too 
small, the slow mode has the dominant viscosity contri- 
bution. At the upper left corner where the layer interac- 
tion is relatively weak, the fast mode is more excited as 
compared to the slow mode. 



In summary, we have calculated the slow relaxation 
rates and the viscoelasticity of a dilute two-layer-vesicle 
solution. We have found the following points: (i) At small 
gap p ss 1, the slowest mode is the squeezing mode. The 
undulation mode appears to be faster, (ii) When the 
inner bilayer radius is small p ks 0, the slow mode be- 
comes the relaxation of the outer bilayer, and the faster 
mode is the relaxation of the inner bilayer. (iii) When 
one decreases p, the relaxation spectrum changes from 
(i) to (ii). (iv) For the complex viscosity, the low fre- 
quency viscosity approaches the hard sphere limit. The 
high frequency viscosity increment is between 12-16 % of 
the hard sphere viscosity increment. The difference be- 
tween the two limits comes from the contributions of the 
two modes, (v) At small gap p > p* , the slow (squeez- 
ing) mode has a small viscosity amplitude, while the fast 
(undulation) mode has a large viscosity amplitude. For 
large gap p < p* , on the other hand, the slow mode has 
the dominant viscosity amplitude, (vi) The crossover size 
ratio p* depends on the interaction between the two bi- 
layers. As Br 2 /n is increased, p* increases from 0.52 
toward unity. 

We have determined the crossover ration p* in Fig. [6] 
The actual bilayer interaction strength B depends on the 
separation d — r 2 (l — p). For the qualitative discussion 



at non-small p, we use the flat layer results in Sec. Ill B 
In Fig. [7] we present both information in the (p, Br 2 j ' n)- 
plot to compare a series of systems with the same outer 
bilayer size r 2 but with different size ratios p. As one 
varies p, the dimcnsionlcss bilayer interaction stre ngth 
Br 2 /n changes according to the DLVO theory Eq. (17) 
(blue lines) or the Helfrich repulsion Eq. ([18]) (red lines). 
If such lines happen to cross the line p* (black line), one 
expects that the two viscosity amplitudes change their 
relative magnitudes. We have indicated such scenarios 
by the dashed lines. 

As shown in the blue dashed line in Fig. [7j an electro- 
statically stabilized system may get large B by having 
a high surface charge density. At small gap, the van 
der Waals attraction may also lower B, causing an inter- 
esting multiple crossing. This means that the viscosity 
amplitudes switch their magnitudes more than once. For 
sterically stabilized system (red lines), we find that the 
dimensionless interaction parameter depends st ron gly on 
the bending rigidity as Br 2 /n ~ k~ 2 (see Eq. (18)). For 
a soft surfactant bilayer of K = 0.75/cbI 1 , we find that the 
slow mode always dominates the viscosity. Whereas for 
a lipid bilayer whose bending rigidity k — 7.5fceT is one 
order of magnitude larger, the red dashed line indicates 
that the squeezing mode is not much excited by shear 
when p > 0.526. 

In this paper, we have only presented the calculation 
for two-layer- vesicles. For MLV with more than two lay- 
ers, the calculation can be performed in the same proce- 
dure, but with a greater algebraic complexity. For MLV 
with N layers, we expect that there are N relaxation 
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FIG. 7. (color online) The dimensionless bilayer interaction 
Br\jK as a function of the size ratio p. The crossover ratio 
p* in Fig. [6] is plotted by the solid black line. The blue lines 
are from the DLVO theory Eq. {TF} with C s = 0.01 M, k d = 
3 x 10~ 8 m, A = 10" 21 J, S = 3 x 10~ 9 m, k B T = 4 x 10~ 21 J. 
The surface potentials are 12 mV and 24 mV for the blue solid 
and dashed lin es, respectively. The red lines are the Helfrich 
repulsion Eq. |l8| with Co = 36/7T 2 , ri — 10" 7 m. The red 
solid and dashed lines are for k = 0.75fc B T = 3 x 10~ 21 J and 
k = 7.5k B T = 3 x 10" 20 J, respectively. 



modes. When the gap is small (r^-i/rN ~ 1), we ex- 
pect that the majority of the N relaxations to bear some 
resemblance to the squeezing mode or the spherical ver- 
sion of the "slip mode" [TT]. For larger gap, on the other 
hand, the relaxations of the N layers may decouple from 
each other. As for the viscoelasticity, we speculate that 
the viscosity amplitudes of the squeezing modes are small 
for weakly interacting MLV. Neutral or weakly charged 
lipid MLV should be an interesting system to investigate 
in this direction. 

In polymer rheology calculation, one may consider a 
step strain for t > 0, where the time relaxation of the 
stress gives the relaxation modulus G(t). This quantity 
can be further converted to the complex modulus G*(u) 
by the Fourier transform. Right after the step strain, one 
often assumes that the polymer deforms in an affine way. 
Do we implicitly use the affine deformation approxima- 
tion for the MLV rheology? A step strain for vesicle solu- 
tion will induce a short but fast bilayer movement, where 
the D terms dominate the left hand side of Eq. (19 1. 



Therefore the initial layer displacements will be propor- 
tional to D _1 • e 2 cx (— Diz/Du, 1) which is compati- 
ble with the incompressible constraints (see Eq. (Bl)) at 
both layers. In the small gap limit, D -1 • e 2 behaves like 
(p 2 , 1). Compared with the affine deformation u n oc r„, 
or in terms of our chosen variables (rfui, r|w 2 ) oc (p 3 , 1), 
it is clear that our theory does not use the affine defor- 
mation approximation. 
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Appendix A: Solution of Stokes equation 

For the incompressible solenoidal flow, the velocity can 
be expressed as 

v = V x (Vt/> x r) + V x (Cr), (Al) 

where the scalar functions ip and £ are the defining func- 
tion for the poloidal and toroidal flow fields, respectively. 
Note that our definition of the defining function differs 
by a factor r compared with the ones in the book by 
Chandrasekhar [25] . 

For an incompressible fluid, the divergence of the pres- 
sure gradient vanishes, i.e., 



V 2 p = 0. 



(A2) 



10 



Therefore the pressure gradient is also a solenoidal field, 
and can be described by the above decomposition. Since 
the toroidal part can be written as (VC) X r, the condi- 
tion dg(d v p) — d v (dgp) requires that dg(— sin# dgC,) = 
c> v (<9 v C/ sin6>) or L 2 ( — 0, where 



Using the standard method of separation of variables, 
the general solution of Eqs. (A2|, (|A8[) and (A9) are 



P = Po + y2(p\„/ +p] 1 m r 1 l )Yi m , 



(A10) 



L = — 



1 o> 
shiflM 



d 



sm t 



i d 2 



d9 J sin 2 9 d V 2 ' 



(A3) 



Hence the pressure gradient is only poloidal, and written 
as 

Vp = Vx (V* x r) 

.i 2 * -1 d d(rV) „ 1 d 9(r#) , A , s 



+ E 

where Yi m (9,tp) are the spherical harmonics. 



n l r l+2 II 

77(2Z + 2)(2Z + 3) 2ry/(2/ - 1) } lmi 1 ' 



r 89 dr 



r sin 9 dip dr 



where \P is the defining function of Vp. Comparing the 
and directions, we can set 



P = 



d(r^) 
dr 



Appendix B: The surface incompressible limit 

The surface incompressibility condition also constrains 
the velocity field as 



(A5) 



2v r (r n ) 



for these two directions. The Laplace equation for the 
pressure Eq. ( A2 1 then implies that 



+ V± • V]_(r n ) = 



(Bl) 



at both r = rj and r 2 . Because the fluid is incompress- 



1 d ( ,<9 2 (r1<) 



r 2 dr 



dr 2 



L 2 9(r*) 
r 2 dr 



(A6) 



ible, i.e., V ■ v = 0, Eq. (Bl ) can also be expressed as the 
equivalent form 



Therefore the radial component of the gradient pressure 
can be either expressed as L 2 ^/r or d 2 (r^>). The latter 



dv r (r n ) 
dr 



= 0, 



(B2) 



expression is consistent with the identification Eq. ( A5 1 



For the velocity field, we will drop the toroidal part 
by setting £ = 0. This is justified because the tan- 
gential force field is a surface gradient and drives only 
the poloidal flow. To obtain the defining function ip, we 
take the curl of Eq. 



at the bilayers. We prefer this condition because it sim- 
plifies the calculation. 

We now consider the relaxation rates of a small per- 
turbation described by the spherical harmonics Yi m (9, ip). 
To simplify the notation, we drop the subscript "Zm" of 
the coefficients ip] m , ipj^, p\ m , and p}^. For the region 



get V x V v = 0. Since within the Mayer i (0 < r < r ^ we set ipf m = and 



V v = — VxVxv holds for incompressible flow, we 
have 



VxVxVxVx[Vx (ipr)} = 0, 



(A7) 



Ira 

Phn = to drop the functions which are singular at r = 0. 
We denote ipj m and p\ m as ip l A and p \ , r espectively. The 
2D incompressibility condition Eq. (B2) relates the two 
remaining coefficients as 



where in the square bracket, an equivalent form of \7ip x r 
is used. As detailed in Ref. [25], each curl will switch 
poloidal and toroidal parts. The double curl will preserve 
the type and modify the defining function by —V 2 . The 
curl of Eq. (15]) becomes V x (rV 4 ip) = or 



2r?(2Z + 3)(Z- 1) 



Pa- 



(B3) 



In terms of the radial velocity amplitude v\, we can fur- 
ther express the pressure coefficient as 



V 4 V> = 0. 



(A8) 



To find the coupling between thepressure p and the ve- 
locity poloidal function ip in Eq. pi, we rewrite — V 2 v = 
VxVxv = Vx [V(-V 2 V>) x r]. Then Eq. (5) is satisfied 
when 

- * + r]V 2 ip = 0. 

Here we prefer to use p instead of vf. Operating the 
above equation by L 2 jr and substituting the combination 
L 2l i>/r by d r p, we obtain 



Pa 



77(2/ + 3)(/-l)r; 
I 



-Vi. 



(B4) 



For the region exterior to the bilayer 2 (r 2 < r < 00), 
the coefficients ip\ m are zero, except ipko, which needs 
to be chosen to give the far flow Eq. ([8]). Comparing 
the radial component v r from Eq. <Js]) and L 2 ip/r, we set 
V4o = r/3. We also set p\ m — so that p does not diverge 
at r — >• 00, and denote ipj^ and pj^ as ipc anc ^ Pc' res P ec ~ 
tively. Then the 2D incompressibility condition Eq. ( B2 ) 
relates the two remaining coefficients as 



dp 
dr 



r]L 2 V 2 ip. 



(A9) 



2t7(2Z - 1)(Z 



-yc 



12 



r 5 2 TS l2 6 m0 . (B5) 
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Using the radial velocity amplitude v 2 , we can express 
the pressure coefficient as 



Pc 



V (2l-l)(l + 2)r l 2 
l + l 



v 2 - 10vr 3 2 TS l2 S m0 . (B6) 



the coefficient of Eq. ( B4 1 with n replaced by r 2 , whereas 



G21 becomes the coefficient of Eq. (B6) (without the T 
term) with r 2 replaced by r%. 

Appendix C: The stress and the tension 
perturbation 



For the region between the two bilayers (r\ < r < r 2 ), 
the expressions are more complex. The incompressibility 



Eq. (B2 1, evaluated at r% and r 2 , provides two conditions 



between the four coefficients 



i> I B = F 11 p I B + F 12 p B l , 
i>B = F 2 ipi + F 22 p I l 



(B7) 



where 



F\2 
i*21 



r 2l+3 _ 21+3 

In I 1 



21 + 1 



2ri(2l + 3){l- 1) (r; 



2r](2l - - 1) 



21+1 21+1 

' 1 ' O 



2r?(2Z + 3)(Z + 2) (r 



2^.2i + l 



i*22 — 



rrr; 



1'2 



22 + 1 
2 

r 2„2/+l 
r 2' 1 



„22+n 



277(2^- l)(Z + 2) (r; 



21+1 
2 



„22+n 



(B8) 



We prefer to use the variables V\ and v 2 instead of p l B 
and Pq. Then the pressure coefficients can be expressed 
as 



Pb = G n vi + G 12 v 2 , 
G 21 v 1 + G 22 v 2 , 



n 11 
Pb 



(B9) 



where 



G11 = ^{(-8l 2 -Al + 12)rf+ 2 r 2 



Gi 



- (4/ + 6)r[ +1 r 2l [l(2l + l)r 2 - (I + 2) (21 - l)r 2 ]} 
77 



/G ( 



{(-8Z 2 -4Z + 12)nr, 



31+2 
2 



(61 + 4)r 2l r l 2 +1 [l(2l + \)r\ - (I + 2)(2l - l)rf] } 



G 2 i = |^ I g{-(2/ 2 + 3/ + l)rf+ 4 r 2 2 ' 

+ (2/ 2 + Z - 3)rf +2 r 2i+2 + (2/ + 4)r 1 +1 r 4 '+ 3 } 



G22 = #^{-(2/ 2 + ^ + l) 



,3Z+4 22 



(/ + 1)G 

+ (2? 2 + z-3) 



,32+2 22+2 



+ (2Z + 4)r 2 + V 



2+1^42+3 1 



with 



G =r 4/ + 4 [Ap + Ap 4l+3 -(2l + l)-p 



-(6-8/-8/ 2 )p 



21+2 



2 p 2l+i 

(21 + l) 2 p 21 } . 



(BIO) 



(BH) 



When the two bilayer are well separated (r x « r 2 ), both 
G11 and G 22 become small. In this limit, Gi 2 becomes 



The radial component of the normal stress appears in 
Eq. ([6]). From its component form a rr — —p-\-2r\d r v r and 
the 2D incompressibility condition Eq. (B2 ), it is just the 



negative pressure. Therefore the stress differences at the 
two bilayers are 

o rr (r\) - a rr (r^) = ^ [(Pa~.PbM -Pb*T' -1 ] Y lm, 

I in 

a rr (r+) - a rr (r 2 ) = ^ [p l B r l 2 + (pg - pg)rj 1_1 ] Y lm . 

lm 

(Cl) 

The large bilayer stretching modulus E suppresses the 
surface density perturbation. For the slow relaxation, we 
will take the limit E — > 00 to eliminate the parameter 
E. Within this limit, both <pi and (f> 2 approach unity, 
so that the tension perturbation 5ji and 5j 2 become the 
Lagrange multipliers. Then the values of the Lagrange 
multipliers are determined by Eqs. ^ and (11) instead 
of their original definitions, as discussed below. 

In Eq. @, we arc interested in the difference between 
the tangential stress across the bilayer. We express the 
tangential stress as a function of ip: 



<j± = 6<r r g + (pa 



Tip 



<9 2 V> 
dr 2 



(L 2 



2 - 
r 



(C2) 



We now replace d 2 ip using Eq. (A9), and limit our dis- 



cussion to I 7^ mode. The velocity field V x (Vip x r) 
has the components similar to Eq. ( A4 1 . We can use the 
r component v r = L 2 ip/r to eliminate %b as (r/L 2 )v r so 
that 



r J2 



a+=V_ 



dp 
J^dr 



2-1 
L 2 



2-q 



r dv r 
E 2 ~dr 



, (C3) 



where L 2 — > 1(1 + 1) for spherical harmonics with nonzero 
I. Because v r is continuous across the bilayer and d r v r 
vanishes on the bilayer, only the first term can be differ- 
ent across the bilayer. This is the only important term 
for the tension perturbation in Eqs. ^ and (11): 



Sjn = - 



We then obtain 



dp 
L 2 



'dp 
L 2 V* 



(C4) 



5 7l = ^(%_^^+i + a rr i)y lro , 



572 = 



1 + 1 

A, 

z + i f2 



r 2l-X + P]tl]& r -l ]Ylnt . 



(C5) 
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